package fem2;

import math2.CrackTipBasis;
import math2.FunctionBasisOnRnToR;
import math2.FunctionRnToR;
import math2.LagrangeFB;
import math2.LagrangeFBOnHex8;
import math2.LagrangeFBOnQ4;
import math2.LagrangeFBOnT3;
import math2.LagrangeFBOnT6;
import math2.LagrangeFBOnTet10;
import math2.LagrangeFBOnTet4;
import math2.PolynomialBasisOnRnToR;
import math2.SerendipityFBOnQ8;
import fem2.enu.MeshPartType;
import fem2.level_set.Crack;

/**
 * singleton class to manage all basis functions for all type of geometry
 * 
 * @author hbui
 * 
 */
public class BasisManager {
	private static BasisManager instance = new BasisManager();

	private BasisManager() {
	}

	/**
	 * Get the instance of BasisManager singleton
	 * 
	 * @return the Singleton instance
	 */
	public static BasisManager getInstance() {
		return instance;
	}

	PolynomialBasisOnRnToR basisT3 = new LagrangeFBOnT3();
	PolynomialBasisOnRnToR basisQ4 = new LagrangeFBOnQ4();
	PolynomialBasisOnRnToR basisT6 = new LagrangeFBOnT6();
	PolynomialBasisOnRnToR basisQ8 = new SerendipityFBOnQ8();
	PolynomialBasisOnRnToR basisTet4 = new LagrangeFBOnTet4();
	PolynomialBasisOnRnToR basisTet10 = new LagrangeFBOnTet10();
	PolynomialBasisOnRnToR basisHex8 = new LagrangeFBOnHex8();
	PolynomialBasisOnRnToR basisLine2 = new LagrangeFB(1);
	PolynomialBasisOnRnToR basisLine3 = new LagrangeFB(2);

	FunctionBasisOnRnToR cb = new CrackTipBasis();

	public FunctionBasisOnRnToR getBasis(MeshPartType mpType) {
		if (mpType == MeshPartType.TRI3) {
			return basisT3;
		} else if (mpType == MeshPartType.QUAD4) {
			return basisQ4;
		} else if (mpType == MeshPartType.TRI6) {
			return basisT6;
		} else if (mpType == MeshPartType.QUAD8) {
			return basisQ8;
		} else if (mpType == MeshPartType.TET4) {
			return basisTet4;
		} else if (mpType == MeshPartType.TET10) {
			return basisTet10;
		} else if (mpType == MeshPartType.HEX8) {
			return basisHex8;
		} else if (mpType == MeshPartType.LINE2) {
			return basisLine2;
		} else if (mpType == MeshPartType.LINE3) {
			return basisLine3;
		} else {
			throw new Error("this mesh part type is not supported");
		}
	}

	/**
	 * @param mp
	 * @param c
	 * @param psi
	 * @param phi
	 * @param r2
	 * @param r
	 * @param dBdr
	 * @param dBdtheta
	 * @param xi
	 * @return
	 */
	private double[] computeCrackTipGradient(MeshPart mp, Crack c, double psi, double phi,
			double r2, double r, double dBdr, double dBdtheta, double... xi) {
		double[] dpsi = c.getLevelSet("psi").gradientAt(mp, xi);
		double[] dphi = c.getLevelSet("phi").gradientAt(mp, xi);

		double[] dr = new double[] { (psi * dpsi[0] + phi * dphi[0]) / r,
				(psi * dpsi[1] + phi * dphi[1]) / r };
		double[] dtheta = new double[] { (dpsi[0] * phi - dphi[0] * psi) / r2,
				(dpsi[1] * phi - dphi[1] * psi) / r2 };

		return new double[] { dBdr * dr[0] + dBdtheta * dtheta[0],
				dBdr * dr[1] + dBdtheta * dtheta[1] };
	}

	/**
	 * create crack tip basis for the mesh part
	 * 
	 * @param mp
	 * @param c
	 * @return
	 */
	public FunctionRnToR[] createCrackTipBasis(final MeshPart mp, final Crack c) {

		FunctionRnToR B1 = new FunctionRnToR() {

			@Override
			public double valueAt(double... xi) {
				double psi = c.getLevelSet("psi").valueAt(mp, xi);
				double phi = c.getLevelSet("phi").valueAt(mp, xi);

				double r = Math.sqrt(Math.pow(psi, 2) + Math.pow(phi, 2));
				double theta = Math.atan2(psi, phi);

				return Math.sqrt(r) * Math.sin(theta / 2);
			}

			@Override
			public double[] gradientAt(double... xi) {
				double psi = c.getLevelSet("psi").valueAt(mp, xi);
				double phi = c.getLevelSet("phi").valueAt(mp, xi);
				double r2 = Math.pow(psi, 2) + Math.pow(phi, 2);
				double r = Math.sqrt(r2);
				double theta = Math.atan2(psi, phi);

				double dBdr = Math.sin(theta / 2) / (2 * Math.sqrt(r));
				double dBdtheta = Math.sqrt(r) * Math.cos(theta / 2) / 2;

				return computeCrackTipGradient(mp, c, psi, phi, r2, r, dBdr, dBdtheta, xi);
			}
		};

		FunctionRnToR B2 = new FunctionRnToR() {

			@Override
			public double valueAt(double... x) {
				double psi = c.getLevelSet("psi").valueAt(mp, x);
				double phi = c.getLevelSet("phi").valueAt(mp, x);

				double r = Math.sqrt(Math.pow(psi, 2) + Math.pow(phi, 2));
				double theta = Math.atan2(psi, phi);

				return Math.sqrt(r) * Math.cos(theta / 2);
			}

			@Override
			public double[] gradientAt(double... xi) {
				double psi = c.getLevelSet("psi").valueAt(mp, xi);
				double phi = c.getLevelSet("phi").valueAt(mp, xi);
				double r2 = Math.pow(psi, 2) + Math.pow(phi, 2);
				double r = Math.sqrt(r2);
				double theta = Math.atan2(psi, phi);

				double dBdr = Math.cos(theta / 2) / (2 * Math.sqrt(r));
				double dBdtheta = -Math.sqrt(r) * Math.sin(theta / 2) / 2;

				return computeCrackTipGradient(mp, c, psi, phi, r2, r, dBdr, dBdtheta, xi);
			}
		};

		FunctionRnToR B3 = new FunctionRnToR() {

			@Override
			public double valueAt(double... x) {
				double psi = c.getLevelSet("psi").valueAt(mp, x);
				double phi = c.getLevelSet("phi").valueAt(mp, x);

				double r = Math.sqrt(Math.pow(psi, 2) + Math.pow(phi, 2));
				double theta = Math.atan2(psi, phi);

				return Math.sqrt(r) * Math.sin(theta) * Math.sin(theta / 2);
			}

			@Override
			public double[] gradientAt(double... xi) {
				double psi = c.getLevelSet("psi").valueAt(mp, xi);
				double phi = c.getLevelSet("phi").valueAt(mp, xi);
				double r2 = Math.pow(psi, 2) + Math.pow(phi, 2);
				double r = Math.sqrt(r2);
				double theta = Math.atan2(psi, phi);

				double dBdr = Math.sin(theta / 2) * Math.sin(theta) / (2 * Math.sqrt(r));
				double dBdtheta = Math.sqrt(r)
						* (Math.cos(theta / 2) * Math.sin(theta) / 2 + Math.sin(theta / 2)
								* Math.cos(theta));

				return computeCrackTipGradient(mp, c, psi, phi, r2, r, dBdr, dBdtheta, xi);
			}
		};

		FunctionRnToR B4 = new FunctionRnToR() {

			@Override
			public double valueAt(double... x) {
				double psi = c.getLevelSet("psi").valueAt(mp, x);
				double phi = c.getLevelSet("phi").valueAt(mp, x);

				double r = Math.sqrt(Math.pow(psi, 2) + Math.pow(phi, 2));
				double theta = Math.atan2(psi, phi);

				return Math.sqrt(r) * Math.sin(theta) * Math.cos(theta / 2);
			}

			@Override
			public double[] gradientAt(double... xi) {
				double psi = c.getLevelSet("psi").valueAt(mp, xi);
				double phi = c.getLevelSet("phi").valueAt(mp, xi);
				double r2 = Math.pow(psi, 2) + Math.pow(phi, 2);
				double r = Math.sqrt(r2);
				double theta = Math.atan2(psi, phi);

				double dBdr = Math.cos(theta / 2) * Math.sin(theta) / (2 * Math.sqrt(r));
				double dBdtheta = Math.sqrt(r)
						* (-Math.sin(theta / 2) * Math.sin(theta) / 2 + Math.cos(theta / 2)
								* Math.cos(theta));

				return computeCrackTipGradient(mp, c, psi, phi, r2, r, dBdr, dBdtheta, xi);
			}

		};

		return new FunctionRnToR[] { B1, B2, B3, B4 };
	}
}
